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ABSTRACT 

Cosmic ray anisotropy has been observed in a wide energy range and at different angular scales by 
a variety of experiments over the past decade. However, no comprehensive or satisfactory explanation 
has been put forth to date. The arrival distribution of cosmic rays at Earth is the convolution of 
the distribution of their sources and of the effects of geometry and properties of the magnetic field 
through which particles propagate. It is generally believed that the anisotropy topology at the largest 
angular scale is adiabatically shaped by diffusion in the structured interstellar magnetic field. On the 
contrary, the medium- and small-scale angular structure could be an effect of nondiffusive propagation 
of cosmic rays in perturbed magnetic fields. In particular, a possible explanation of the observed 
small-scale anisotropy observed at TeV energy scale, may come from the effect of particle scattering 
in turbulent magnetized plasmas. We perform numerical integration of test particle trajectories in 
low-/3 compressible magnetohydro dynamic turbulence to study how the cosmic rays arrival direction 
distribution is perturbed when they stream along the local turbulent magnetic field. We utilize 
Liouville’s theorem for obtaining the anisotropy at Earth and provide the theoretical framework for 
the application of the theorem in the specific case of cosmic ray arrival distribution. In this work, we 
discuss the effects on the anisotropy arising from propagation in this inhomogeneous and turbulent 
interstellar magnetic field. 

Subject headings: cosmic rays - magnetic fields - magnetohydrodynamics (MHD) - scattering - turbu¬ 
lence 


1. INTRODUCTION 

Cosmic rays are found to possess a small but measur¬ 
able anisotropy in their arrival direction distribution at 
Earth. The origin of the observed anisotropy is not yet 
understood. However, it is reasonable to assume that 
it is a combination of effects correlated to the distribu¬ 
tion of the galactic sources of cosmic rays, the geometry 
and turbulence properties of the galactic magnetic field, 
and the propagation in interstellar magnetized plasmas. 
These are likely to be responsi ble for the complex shape 
of the energy spectrum as well (Gaisser et al. 2013). Since 
we don’t know the locations of cosmic ray sources and 
we lack details of the interstellar magnetic field, under¬ 
standing these observations of anisotropy is not an easy 
task. 

Above the energy range where cosmic rays are directly 


affected by inner heliospheric processes ( see, e.g., Elorin- 


ski et al. (|2013|); |Manuel et al.| (|2014|); jk'lor inski et al 
(|2015|)), a statistically significant anisotropy has been 
observed by a variety of experiments, sensitive to dif¬ 
ferent energy ranges (from tens of GeV to a few PeV), 
located on or below the Earth’s surface in the North- 


ern Hemisphere (|Nagashima et al.|1988l Hall et al. 

19991 

Amenomori et al.||2UU5l |2UU6| |Guillian et al. 

WTT 

Abdo 

et al. 112009 |Aglietta et al.||2009| Zhang et al 

.1120091 Mu- 

nakata et al.| 20101 lAmenomori et al.||20111 

de Jong et 

al. 2011 Shuwang et al. 

20111 Bartoli et al, 

2015 

) and 

in the Southern Hemispi 

lere (lAbbasi et al. 2010, 

2011, 

2012 Aartsen et al.||2013 

). 


trivial fashion. Erom about 100 GeV to tens of TeV, 
it has been observed to have an approximately consis¬ 
tent structure at the largest scale, although its ampli¬ 
tude appears to increase with energy. Above a few tens 
of TeV, the observed progressive change in the anisotropy 
topology may indicate a transition between two processes 
shaping the particles’ arrival distribution at Earth. The 
observation could be qualitatively explained on the ba¬ 
sis of diffusive propagation of cosmic rays in the Milky 
Way from stochastically distributed sources, responsible 
for generating a gradient in cosmic ray density. 

Numerical studies of particle propagation in a scenario 
of homogeneous and isotropic diffusion in the Galaxy 
predicts that the cosmic ray density gradient, and the 
consequent induced anisotropy, has a dipole shape with 
direction toward the source of the particle and with am¬ 
plitude that directly depends on the diffusion coefficient. 
In particular, for a given realization of galactic source 
spatial distribution, the dipole anisotropy wou ld point 
toward the source with the largest contrib ution (|Erlykin 


fc: Wolfendale||2QQ6| Blasi fc Amato||2Q12| |Ftnskin||2UI^ 

Pohl fc Eichler||2QT3 Sveshnikova et al.||2Q13 Savchenl^ 

et al.|2015 ), which may change with energy, in agreement 

with observations. On the other hand, the systematic 
overestimation of the anisotropy amplitude may be par¬ 
tially compensated by the fa ct that diffusion is expecte d 


Effenberger et al. ( 2Q12| )) 


to be anisotropic (see, e.g. 
thus modifying the expected cosmic ray density' gradient 
shape as a function of the source dire ction with respect 
to th e regular galactic magnetic field (Kumar & Eichler 
2014). The misalignment between the cosmic ray density 
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gradient and the regular galactic magnetic field would 
prevent pointing to any specific source, and it would sup¬ 
press the anisotrop y amplitude to a value closer to what 
has been observed (Mertsch & Funk |2Q15 ). 

Another scenario is that it is the transition from 
heliospheric- to interstellar-dominated contributions, 
starting at a 10 TeV energy scale and culminating around 
1Q Q-2QQ TeV, at the or i gin o f the shift in anisotropy. 
In Desiati & Lazarian (2013), this scenario was pro- 
posed noting that 10 lev protons have an average gy- 
roradius, in a /iG scale magnetic field, on the order 


of the thickn ess of the heliosphere (see, e.g., Pogorelov 
et ah ([2006)). In addition, the dynamical instabilities 
of the heliospheric magnetized plasma at smaller scales 
may generate strong scattering that redistributes the ar¬ 
rival direction of TeV cosmic rays. This heliospheric 
sce nario is studied and presente d in our companion pa¬ 
per Lopez-Barquero et al. (2016). Such strong scattering 


may be able to produce large localized particle gradl 
ents, experimentall y interpreted as medium - or small- 
scale anisotropy. In Schwadron et al. (2014), a scenario 
of weak influence of the heliospheric magnetic field on 
TeV protons was explored, thus interpreting the obser¬ 
vations as directly related to the ordering of the loc al 
interstellar magnetic field (see also Zhang et al. (2014|)). 

The anisotropy appears to possess a complex angular 
structure with evidence of a harder cosmic ray spectrum 
within the localized excess region in the apparent direc¬ 
tion of the heliospheric tai l (Am nnmnorimtah [2007 jAbdo 


2008 


Bart oh et al 


2013[ Abeysekara et al. 


TUTA) 


4'he deconiposition of 'I'eV cosmic ray anisotropy in the 
individual spherical harmonic contributions shows that 
the arrival direction distribution is dominated by large- 
scale structures (such as, e.g., dipole and quadrupole) 
with a relative intensity on the order of 10“^, but that 
medium- and small-scale angular structures are signifi¬ 
cantly contributing with relative intensities below 10“^. 

The experimental determination of small angular scale 
anisotropy is generally performed by filtering out the 
large-scale modulations from the observed arrival direc¬ 
tion distribution, thus retaining all structures with large 
angular gradients. Some of the small-scale anisotropy 
features seem to be correlated to regions in the sky where 


the globa l aniso tropy has large variations (see Desiati & 
(2013)) or may be an effect of re-acceleration 


Lazarian 


heliosphere ( 

Lazarian & Desiati[ 

[2010 Desiati & Lazar- 

ian 

2012 

). However, globally, tJ 

tie observed small scale 


and, therefore, possibly a natural consequence of cosmic 
ray propagation in the local turbu lent magnetic field i n 


the presence of a global anisotropy (Giacinti & Sigl 2012). 
The global anisotropy, at all angular scales, arises from 
the same physical processes, thus it is impossible to dis¬ 
entangle its origin. However, as a first approximation, it 
is generally assumed that the anisotropy at the largest 
scale is dominated by global physical processes (such as a 
density gradient from sources of cosmic rays or from con¬ 
vective effects originated by large-scale cumulative stellar 
winds, and from propagation through the regular galac¬ 
tic magnetic field), while the small angular scales are 
dominated by local processes. 

A complex angular power spectrum is asymptotically 


generated by progressive decomposition of the energy of 
an initial anisotropy distribution (for instance a dipole) 
into higher multi poles, by the effect of scattering off mag - 
netic turbulence (|Ahlers||2014t [Ahlers & Mertsch||2015|). 
The conservation of phase space density, as stated by 
Liouville’s theorem, predicts, for an idealized situation 
of a homogenous large-scale anisotropy, the total sum of 
the multipoles’ angular terms is conserved. This makes 
it possible to ge nerate small-scale structures, as shown 


m 


Ahlers (2014). 

Turbulence m astrophysical plasmas have significant 
effects on particle propagation, in that the stochastic na¬ 
ture of magnetic field lines is transmitted to the parti¬ 
cles’ trajectories. If particles are tied to magnetic field 
lines, the maximum perpendicular diffusion rate is set 
by the rate of perpendicular field line wandering, scaled 
by particle velocity (field line ra ndom walk) . This has 


been extensively d iscussed by, e.g. Jokipii (|19 66|);|J o kipii 
& Parker! (|1969|); [Rechester & Rosenbluthj (|1978D; 

- M - WISE 


acalone & Jokipii ( |1999| ); Minnie et al. ( |2UQ* 




cEi 


in the context of dinusion at distance scale larger 
than the turbulence injection scale. On scales smaller 
than the injection scale, particles are characterized by 
super-diffusion in the perpendicular direction of the 
mean magnetic field. This behavior is described by the 
so-called Richards on diffusion, where par ticl e separation 


TOWS as (time)^/^ ( [Lazarian fc Yan 2014). InEyinketal 


2011[) it was shown that this is directly connected to the 
separation of tu rbulent magnetic field lines , which grows 
as (distance)^/^ ([Lazarian & Vishniac|l999). In this case, 
the stochastic nature of the astrophysical magnetic fields 
(such as, e.g. the interstellar magnetic field, or, at larger 
scales, the intercluster magnetic field), inevitably pro¬ 
duces chaotic particle trajectories, meaning that the ge¬ 
ometry of particle trajectories is highly sensitive to the 
actual initial conditions. Diffusion by magnetic field line 
wandering is found to be a dominant contribution and 
stronger than the extreme case of Bohm diffusion, where 
scattering frequency is one per particle gyration. The im¬ 
portant aspect that determines the properties of a large 
ensemble of particles is its statistical nature, which in this 
case is influenced by the properties of the magnetic field 
and specifically by the induced scattering rate. While 
individual trajectories may have chaotic properties and 
are, therefore, practically/realistically unpredictable, the 
large ensemble of particles can still be deterministically 
described. For these reasons, the recorded spatial distri¬ 
bution of the ensemble will be statistically determined 
and a direct consequence of the properties of the turbu¬ 
lent magnetic field. In this work, a study of particle prop¬ 
agation in compressible magnetohydrodynamics (MHD) 
turbulence is performed in the context of its effects on 
arrival direction distribution. 

In addition, depending on the degree to which mag¬ 
netic field lines diverge on scales less than the particle 
gyroradius, pitch angle scattering on small-scale mag¬ 
netic perturbations affects particle distribution at the 
large spatial scale, thus increasing the diffusion coeffi¬ 
cient depending on the large geom etrical scale of the 
magne tic field turbulence (see, e.g., Desiati & Zweibel 
(2014) and references therein). The anisotropic nature 
of interstellar turbulence and the properties of turbu¬ 
lence itself can significantly complicate the description 
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of cosmic ray transport (see, e .^., |Yan fc Lazarian| ( |2QQ8 
'2011); Beresnyak et al. (2011|)). 

Tnis paper is organized as follows. In section we 
describe the turbulent magnetic field used and the test 
particle trajectory integration. In section we discuss 
the validity of applying Liouville’s theorem in the context 
of this work by statistically assessing the level of conser¬ 
vation of magnetic moment. The results of the study are 
presented in section]^ and discussed in sectionCon¬ 
cluding remarks follow in section 

2. COSMIC RAY PROPAGATION IN TURBULENT 
MAGNETIC FIELDS 

Cosmic rays, which are accelerated in their sources 
and “injected” into the interstellar medium (ISM), are 
free to propagate through the interstellar magnetic field. 
Globally, the galactic magnetic field is characterized by a 
large-scale regular co mponent and a sniall-scale r andom 
component (see, e.g., |Jansson & Farrar | (|2012a|b|)). The 
regular component can be described as a superposition 
of spiral and toroidal structures, possibly with a contri¬ 
bution perpendicular to the galactic plane, and the ran¬ 
dom component represents the stochastic perturbations 
of the regular field caused by the dynamics of the galaxy 
and its density distribution. An important property of 
the interstellar magnetic field is turbulence. A variety of 
observations show that the coherence scale of turbulent 
magnetic fields is on the order of 10 pc in spiral arm re¬ 
gions (with more frequent stellar formation ac tivity) and 
on the orde r of 100 pc in the interarm regions (Haverkorn 
et al.||2008|). The injection scale of turbulence is deter¬ 
mined by the scale at which the magnetic perturbation 
is generated (for instance, by stellar collapse or binary 
mergers). 

Astrophysical plasmas are typically highly ionized and 
have high Reynolds numbers, thus the dynamics of the 
flow is dominated by nonlinear convective processes at 
the largest scale. In such conditions, turbulence de¬ 
velops and magnetized Alfvenic eddies dynamically cas¬ 
cade to smaller scales and progressively elongate along 


the magnetic fie l d line s , as initially proposed by Srid- 
har fc Goldreich|(|'l9 94|);|Goldreich fc Sridhar](11995) (see 
also|Lazarian|(|2007|) for a review). T'his mod^ of incom¬ 
pressible IVIHL) turbulence predicts a Kolmogorov-type 

energy power spectrum oc in terms of the 

wave-vector component perpendicular to the local direc¬ 
tion of the magnetic field, while the parallel component 

of the wave vector is k\\ = k^/^. In those pioneering 
papers, the theory assumes the injection of energy at 
scale L and the injection velocity equal to the Alfven 
velocity in the fluid Va, i.e., the Alfven Mach number 
Ma = {VlIVa) = 1 (i.e., trans-Alfvenic turbulence), 
where Vl is the plasma velocity at injection scale L. The 
model was later generalized for both sub-Alfven ic (i.e., 
Ma < 1) and super-Alfvenic (i.e., Ma > 1) cas es (iLazar- 
ian fc Vishn iac 1999| Lazari an |2QQ 6 ) (see also Lazarlan 
et al. ( |2U12[ ); [Burkhart et aL| ([^14|) ). Typically, the IblVl 


is characterized by Ma ^ 1 ([Gaensler et al. [2011). This 
means that magnetic field lines do not typically fluctuate 
too far from the mean direction. 


Alfven, slow, and fast waves. Alfven modes are incom¬ 
pressible, while slow and fast modes are compressible. 
The compressible modes are conjectured to resemble in¬ 
compressible behavio r at high-/3 (i.e., high gas to mag¬ 
netic pressure ratio) (jLithwick & Goldreich||2001 ); how- 
ever, ISM plasmas are typically characterized by low-;5. 
Gho & Lazarianj (|2002|) investigated the scaling proper- 
ties of low-/i couipressible sub-Alfvenic MHD turbulence 
and confirmed that slow compressible modes behave as 
the Alfven incompressible modes but also that the fast 
modes are isotropic, since their velocity does not depend 
on magnetic field direction. 

The spatial distribution of particles propagating in a 
turbulent field is affected by the magnetic perturbations 
within the particle autocorrelation length scale (or mean 
free path). To study the correlation between turbulence 
and cosmic ray distribution, trajectories of test parti¬ 
cles were integrated in an MHD turbulent magnetic field. 
The methodology used is described in the next two sec¬ 
tions. 

2.1. Turbulent magnetic field 

Test particle trajectories are integrated in a compress¬ 
ible sub-Alfvenic isothermal MHD turbulence i n low-/3 
based on nume rical calculations developed by [Gho fc 
Lazarian (2002). In the model, turbulence is driven 
solenoidally m Fourier space, setting the velocity and 
density fields initially to unity. The calculation is per¬ 
formed in a cube with side Li)ox = 512 grid-points, 
with inertial range from Li^j = 204.8 grid-points (i.e., 
0.4 X Lijox) down to a “damping” scale of Lr^in = 5 grid- 
points. 

The average magnetic field is directed along the x- 
axis and turbulence is characterized by a gas-to-magnetic 
pressure value of p ^ 0.2 and by Alfvenic Mach number 
Ma = 0.773. Such a Mach number corresponds to the 
fluctuations being on the order of the mean magnetic field 
at the injection scale, which is in agreement with what 
we would have expected from the local ISM. The exter¬ 
nal mean magnetic field is the only controlled parameter 
in this MHD model. 

2.2. Cosmic ray propagation 

The analysis is performed by integrating proton tra¬ 
jectories in a static magnetic field, using the set of 6- 
dimensional ordinary differential equations 


II 

X 

( 1 ) 

d.r 


( 2 ) 


describing the Lorentz force and the particle velocity u, 
where r is the particle position vector and p the mo- 
memtum. For H, we use one steady st ate r ealization of 
the magnetic field described in section 2.1 The equa- 


While the model by Goldreich & Sridhar (1995) de¬ 
scribes incompressible MHD turbulence, modeling com¬ 
pressible turbulence turned out to be more complex. In 
isothermal plasmas, there are three types of MHD waves: 


tions are integrated using the Bulirsch-Stoer integration 
method with adaptive time step. At each integration 
step, the magnetic field is interpolated using a 3D cubic 
spline, and integration is stopped when particles cross 
the border of the MHD simulation box. The accuracy 
of the particle trajectory integration for this study was 
assessed and is discussed in Appendix The choice of 
one specific realization of the magnet icneld is justified 
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TABLE 1 

Physics parameters of simulation sets 


Set 

Ep 

tl 

^inj 

^mfp 

Particles 

(Bo) 

I 

750 TeV 

0.24 pc 

10 pc 

5.0 pc 

3 X 10® 

3 i^G 

2 

7.5 PeV 

2.4 pc 

100 pc 

50 pc 

10® 

3 /rG 

3 

30 PeV 

10 pc 

100 pc 

60 pc 

10® 

3 fiG 


by the fact that particle velocity is much larger than the 
plasma Alfven velocity, thus induced electric fields can 
be neglected. 

The magnetic field can be interpreted as a snapshot 
of the local interstellar magnetic field. Since the MHD 
simulation can be scaled to an arbitrary injection scale, 
in this study, two physical scenarios are investigated: an 
injection scale of = 10 pc, which is approximately 
the magnitude typical of the Milky Way spiral arms, and 
Linj = 100 pc, which is characteristic of the interarm re¬ 
gions (see section]^. Although the ISM surrou nding the 
sol ar system se ems to have interarm properties ( [Frisch et 


ar]|2012[ 2014), both cases are taken into consideration. 


Xs^summg a particle gyroradius of 5 grid-points (corre¬ 
sponding to the smallest spatial scale in the MHD simu¬ 
lation), the injection scale sets the particle energy scale. 
In order to avoid grossly underestimating the effects of 
magnetic perturbations smaller than the particle gyrora¬ 
dius, the latter is set to 20 grid-points as well. Table 
shows the three trajectory integration data sets used in 
this study. They cover an energy range from 750 TeV to 
30 PeV, partially overlapping the recent cosmic ray ob¬ 


servations reported by the IceCube O bservatory (Abbasi 
et al.||2012| [Aartsen et al.||2013[ |2015|) . 

In order to study the ettect of interstellar magnetic 
turbulence on the arrival direction distribution of cos¬ 
mic rays, a large number of particles would need to be 
injected with randomly uniform directions on a spheric 
surface centered at the Earth and with a radius larger 
than the mean free path. Unfortunately, such method 
would be highly inefficient because a large fraction of 
the injected particles would never reach Earth. The al¬ 
ternative approach is typically to use the so-called “back- 
propagation method,” where the trajectories of particles 
with opposite charge are integrated from Earth, with 
initial directions uniformly distributed, backward into 
outer space. Since energy losses are negligible for proton 
particles, their energy is conserved, and therefore their 
trajectories can be time-reversed, provided there are no 
collisional scattering processes and no resonant scatter¬ 
ing. Under such general conditions, Liouville’s theorem 
states that particle density in phase space is conserved 
along particle trajectories. In the presence of turbulence, 
magnetic fields can vary faster in space than the particle 
gyroradius, thus breaking adiabaticity and effectively in¬ 
ducing collisional processes. In section the applicabil¬ 
ity of Liouville’s theorem in the present case is discussed. 

Eor set I of Table 0 3 X 10^ proton trajectories are in¬ 
tegrated, while for sets 2 and 3, the number of particles 
is 10^. The numerical trajectory integration starts from 
a point at the center of a 3x3 MHD simulation box, with 
uniform direction distribution, and stops when particles 
reach the edge of the integration volume. The integration 
box corresponds to a distance scale of 75 pc (for set I in 
table and of 750 pc (for sets 2 and 3 in Table [^. The 
back-propagated trajectories calculated in this way pro¬ 
vide the information on the effects of the magnetic field 


on an initially uniform particle distribution emanating 
from one point. In other words, they provide a map of 
regions that are magnetically connected to the origin at 
a given distance. The calculation implicitly takes into 
account the dynamical processes of a particle’s motion 
in a magnetic field, including drifts and pitch angle scat¬ 
tering. Under the hypothesis that such trajectories are 
time-reversible (see section]^, they can be interpreted as 
directly propagating from outer space back to the origi¬ 
nal point (assumed to coincide with Earth). This means 
that the particle distribution far from Earth, resulting 
from the back-propagation numerical calculations, cor¬ 
responds to a perfectly uniform arrival distribution at 
Earth. 

The numerically calculated trajectories can be used, 
therefore, to determine the arrival distribution at Earth 
as a consequence of a global anisotropy at a large dis¬ 
tance. Such a global anisotropy, which is the effect of 
a particle density gradient induced by sources of cosmic 
rays or by convective processes, is treated as a weight 
on the forward-inverted trajectories, and it effectively 
produces a nonuniform arrival direction distribution at 


Earth, which is described in section 4.2 


3. THE VALIDITY OF LIOUVILLE’S THEOREM 

Generating a large number of particle trajectories that 
pass through a point in space is implicitly highly ineffi¬ 
cient, since most particles will bypass the target point. 
One possibility is to increase the size of the target to 
record those particles that pass “nearby,” or to appeal 
to Liouville’s theorem. 

Liouville’s theorem states that the density of particles 
in the neighborhood of a given system in phase s pace is 
constant if restric tions are imposed on the system (Gold- 
stein et al.||2002|), as shown below. To determine its va- 
lidity for our specific case of cosmic ray arrival distribu¬ 
tion, w e shall start wit h the continuity equation in phase 
space (Bradt 2008a|b), under the assumption that the 
number of particles stays fixed: 


% = -^- - Vp • ^pp) 


(3) 


= _pV ■{v)-v- V{p) - pVj, ■{F)-F- Vpip) (4) 


where Vp is the del operator in momentum space, p the 

distribution function, and F the applied external force. 

Upon expansion of this expression, we arrive at: 

dt 

The first term on the right-hand side of Eq. van¬ 
ishes, since the divergence of the velocity in phase space 
is zero. The third term in the right-hand side of Eq. 
with Vp • (F ), gives us restrictions on the forces that can 
be applied to the particles. Eor this term to go to zero, 
we need the forces to be conservative and differentiable. 
In particular, they must be p-divergence free. Evidently, 
if collisions are present, they will not fulfill these require¬ 
ments. Now, we can analyze the case of magnetic forces 
and can express this term explicitly: 

Vp ■ {F) = qB ■ (Vp xv)-qv- {Vp x B) (5) 

Assuming that the magnetic field is independent of the 
velocity of the particle, the p-curl of B goes to zero. 
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In this case, the assumption is valid since the particles 
are moving so fast that the possibility of changing the 
magnetic field is negligible. For the p-curl of F, it just 
vanishes, since the velocity is the gradient of a scalar 
function, the relativistic energy. Therefore, if we have 
only pure magnetic forces, the Vp • (F) cancels. 

We arrive at the equation: 


dp 

dt 


+ • V (p) + F • Vp(p) — 0 


(6) 


But this is just the expression for the total derivative 
of the distribution. Therefore, we can reexpress it as: 


? = 0 

dt 


( 7 ) 


which is the precisely the expression for L iouville’s theo¬ 
rem ( Goldstein et al.||2QQ2 


pression tor L v 
(Bradt|[2QQ^ . 


This derivation is for a pure magnetic force, but in fact 
when calculating the particles’ trajectories in a turbulent 
magnetic field, a variety of factors come into play. The 
most significant effect is when particles encounter a re¬ 
gion where the magnetic field varies abruptly, i.e., the 
scale of variation of the magnetic field is shorter than 
the gyroradius of the particle. In this scenario, the tra¬ 
jectory does not have time to adjust smoothly to this 
change, and the interaction can be effectively considered 
a collision. For such cases, the right-hand side of Eq. 


which 


can be modified by the addition of a term, 

takes into account collisions of various origins that are 
differentiated by their exact functional form, given the 
fact that they will cause a n onzero time rate of chang e 


in the distribution function (Baumjohann et al. 1996|). 
Therefore, under these conditions, Liouville's tneorem 
can’t be applied. To ensure that the abruptness in the 
trajectory is limited, we can calculate how adiabatic this 
change is. Having established this, it will ensure that a 
smooth variation will not greatly modify the density in 
phase space. 

In the presence of collisions, the magnetic moment of 
the gyrating particles changes. Therefore, to check for 
the adiabaticity of the trajectories, we can calculate the 
magnetic moment for each particle at each time step and 
find out if it statistically behaves as an adiabatic invari¬ 
ant. 

The relativistic magnetic moment (also called first adi¬ 
abatic invariant) is given by: 


= 


2m\B\ 


(8) 


where p± is the momentum perpendicular to the mag¬ 
netic field B and m the particle mass. This quantity, 
relating magnetic field and perpendicular momentum of 
the particle, is conserved if the field gradients are small 
within distances comparable to the particle gyroradius. 
Using the integrated particle trajectories from the data 
sets of Table the magnetic moment Eq. was calcu¬ 
lated at each integration time step and histogrammed in 
time step slices. In each time slice, the mean value p of 
the magnetic moment of the particle ensemble from each 
data set and the corresponding standard deviation 
were calculated. Eigure shows the evolution of p over 


the integration time of the 30 PeV set (on the left) of Ta- 
blef^and the ratio cr^/p for the three sets (on the right) of 
Tamely A perfect distribution with cr^/p = 0, indicates 
that the particles most likely stay in the same magnetic 
field line, which is unrealistic for a particle moving along 
a turbulent magnetic field. Nonetheless, a distribution 
peaked at a value much larger than one will tell us that 
the particles suffered strong variations in their trajecto¬ 
ries and collision-like interactions happened. The dis¬ 
tributions calculated for our three different energy data 
sets peak at around 0.5 and do not appear to have trends 
or large variations during integration time (compared to 
cr^), as shown on the left of Eigure This indicates 
that even though the particles have interacted with the 
turbulent field, and it has changed their trajectories, the 
changes are relatively smooth and with limited statistical 
impact on the overall particle ensemble. 

Note that the width of the distributions in a^/p means 
that at some level the magnetic moment of the particles 
fluctuates during propagation. Different particles follow 
different magnetic field lines and experience different in¬ 
teractions. There might be also a contribution from the 
limited accuracy of the numerical integration. However, 
as discussed and shown in appendix in the present 
application, the effects due to accuracy limitations are 
not sufficiently large to significantly violate adiabaticity. 
The highest level of possible inaccuracy, where particles 
are stochastically re-distributed at each gyration, pro¬ 
vides a numerical diffusion at the level of Bohm diffusion, 
where particles are scattered every gyration. It has been 
proven that in astrophysical turbulence Bohm diffusion 
is always much smaller t han diffusion induced by mag¬ 
netic field line wandering (jLazarian & Yan|2014|) . So the 
fact that the accuracy of the used trajectory integration 
method is significantly below the Bohm diffusion level 
(see appendix]^ poses no problem on the statistical re¬ 
sults obtained m this study. 

Based on the fact that the ensemble-average {a^/p) < 
1, especially that only magnetic forces are explicitly 
taken into account and that the changes in the parti¬ 
cles’ trajectories are relatively smooth, it is assumed that 
Liouville’s theorem is applicable (see further discussions 
in Section]^. So in this study, back-propagation of the 
particle trajectories is justified in this statistical sense. 
The situation changes if scattering rate is higher, such as 
in the heliospheric mag netic field case studied in [Lopez-"] 
Barquero et al. (2016), where resonant scattering pro- 
cesses produce larger deviations from adiabaticity. 

4. RESULTS 

This section shows the results obtained with t he n u- 
merical calculation data sets described in Section 

4.1. Mean Free Path 

Using the particle trajectory data sets from Tableit 
is possible to evaluate general properties that the ensem¬ 
ble of particles have after a sufficiently long duration of 
propagating in the turbulent magnetic field. The mean 
free path Xmfp is the distance at which the instantaneous 
particle directions become uncorrelated with respect to 
those at time t=0. At this distance, and associated to 
scattering time scale, particles have lost memory of the 
direction distribution at initial conditions. This is a cu¬ 
mulative property of all particles, and it can be estimated 
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Fig. 1.— Left: The mean magnetic moment J1 for the data set at 30 PeV (red circles) with the corresponding moving average with subset 
of 30 time steps (blue line). The magnetic moment is calculated for all the particles at each time step, and then the average for the total 
set is calculated at each step. Right: Histogram of standard deviation of magnetic moment (7^ over mean magnetic moment fL for the three 
data sets of Table ^ The red, blue, and green lines represent 750 TeV, 7.5 PeV, and 30 PeV protons, respectively. The magnetic moment 
is calculated for each particle at all time steps. The mean value and the standard deviation are for the total trajectory. Note that the 750 
TeV set has three times the number of particles as the other sets, as described in Table 
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Fig. 2.— Sky map of arrival direction distribution of the 7.5 PeV proton set of Table Rafter propagation for a distance of 40 pc. The 
decomposition of the initial dipole distribution is shown. On the left is the sky map obtained after time inversion, and on the right is the 
same map after subtracting the dipole component from the map on the left. A Gaussian smoothing with a = 3° was used. 


by calculating the mean distance at which the direction of 
each particle has an angle of 90° from its initial position. 
This definition is equivalent to evaluating the velocity 
autocorrelation function and estimating the distance at 
which it is sufficiently close to zero (i.e., correlation to 
initial condition is lost). 

The results of Xmfp calculations are given in Table 
As expected in this regime, the mean free path increases 
with energy. For the interarm region, where the injec¬ 
tion scale is on the order of 100 pc, Xmfp ^ bO pc for 
30 PeV protons. For an energy four times smaller, 7.5 
PeV, the mean free path decreases to a value of 50 pc. 
In the spiral arm, with injection scale on the order 10 
pc, our calculation for 750 TeV protons is 5.0 pc. Note 
that in the two lowest energy data sets considered here, 
the particle gyroradius is on the order of the damping 
scale of the turbulent field. Therefore, pitch angle scat¬ 
tering arising from smaller scale magnetic perturbations 
is significantly reduced, and may result in the mean free 
path being overestimated to some degree. In the present 


work, the intent is to show the effects of turbulence in 
the particle spatial distribution in relation to the mean 
free path; therefore, whether the value is strictly correct 
is of limited importance in this context. In future stud¬ 
ies, the need to use MHD simulations in a wider inertial 
range will be considered in more detail. 


4.2. Sky Maps of Arrival Direction Distribution 

The particle trajectories numerically integrated with 
Eqs. 1-2 for the sets in Table are used to study the 
properties of their arrival direction distribution that re¬ 
sults from the scattering off magnetic turbulence from 


a particle density gradient. As described in Section [2^ 
the procedure used for determining the sky maps of the 
particles’ arrival distribution makes use of the trajectory 
integration backward in time, uniformly emanating from 
one point assumed to be Earth, until they exit the in¬ 
tegration box. At a sufficiently long distance from the 
origin, particle trajectories accumulate in the direction 
of the mean magnetic field, since the perpendicular diffu- 
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Fig. 3. — Sky maps of arrival direction distributions of 30 PeV 
protons in equatorial coordinates, with the dipole density gradient 
weight at different distances: = 10 pc, 20 pc, 60 pc, and 90 pc 

(from top to bottom). Gaussian smoothing with a = 3° was used. 
On each map, a dipole fit was performed and the resulting dipole 
component was subtracted. 

sion is smaller than that which is parallel to the magnetic 
field lines. A sphere of radius R with center at the origin 
of the back-propagated trajectories is used to record the 


position and velocity direction of each trajectory. The 
trajectories are inverted in time from those locations on 
the sphere to the origin, by virtue of Liouville’s theorem 
(see Sectionj^for a discussion of the validity of Liouville’s 
Theorem ana the consequences of its use in the context 
of the present study). 

It is likely that there are several sources of cosmic 
rays in the Milky Way, perhaps from different popula¬ 
tions injecting particles into the ISM with different en¬ 
ergy spectral shapes. It is also possible to suppose that 
a single source may dominate the cosmic ray distribu¬ 
tion at Earth, depending on the energy range (see Sec¬ 
tion f^. In a scenario of isotropic diffusion, the cosmic 
ray density gradient is expected to be described by a 
simple dipole distribution. This is a similar distribution 
as would be expected if convective transport were the 
dominant source of the cosmic ray density gradient, su ch 
as in the case of superbubbles (Biermann et al. |2Q13|) or 
the effect o f Loop I shell in the local ISM ( [Schwadron 


et al.l 2014). In the general and more likely scenario of 


anisotropic diffusion, cosmic ray propagation along the 
magnetic field line is faster than the perpendicular, thus 
producing a cosmic ray gradi ent ordered along the reg ¬ 
ular magnetic field (see , e.g., Kumar & Eichler (2014); 


Mertsch & Eunk (2015)). Although the density gradient 
is not expected to be a dipole but rather a distribution 
that depends on the ratio of perpendicular to parallel dif¬ 
fusion, in this work it is assumed a simple dipole distri¬ 
bution, regardless of the origin of the underlying density 
gradient of the cosmic rays. 

If the underlying distribution of cosmic rays is perfectly 
uniform, the effects of scattering off magnetic turbulence 
shuffles one isotropic distribution into another isotropic 
distribution. However, with an existing particle density 
gradient, the scattering processes redistribute particles 
from the region of higher density to that of lower density, 
and vice versa, thus creating a complex arrival distribu¬ 
tion that can be described with higher order multipole 
terms of the spherical harmonic expansion. 

In the process of inverting time on the numerically cal¬ 
culated trajectories, a dipole gradient distribution, as¬ 
sumed to be aligned with the direction of the mean mag¬ 
netic field of the MHD snapshot, is introduced as a weight 
on each trajectory at the crossing point on a sphere of 
radius R. The weight is calculated using the angle of 
the particle direction at the crossing point from that of 
the density gradient. The arrival distribution of these 
forward-propagated trajectories at Earth (i.e., the ori¬ 
gin) is then determined. One key factor to remember is 
that only the small-scale angular anisotropy is studied, 
since this is the one that arises from the specific inter¬ 
action with the turbulent magnetic field. Therefore, the 
large scale component, specifically the dipolar compo¬ 
nent of the map at Earth, is fitted and removed. Such a 
dipoie component can have a different direction and am- 
phtude than that injected. Eigurej^ shows the sky map, 
in equatoriai coordinates, of the arrivai distribution of 
the 7.5 PeV protons at Earth before (on the ieft) and 
after (on the right) dipoie subtraction, thus emphasiz¬ 
ing the smaii scaie features. A Gaussian smoothing with 
cr = 3° was used. The residuai map shows medium- and 
smaii-scaie anguiar structures arising from the breakout 
of the underiying dipoie anisotropy after a propagation 
of i? = 50 pc. 
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Fig. 4.— Sky maps of arrival direction distributions of 750 TeV protons (on the left) and 30 PeV (on the right) in equatorial coordinates 
and at propagation distance corresponding to the mean free path. Gaussian smoothing with a = 3° was used. On each map, a dipole fit 
was performed and the resulting dipole component was subtracted. 


Th e maps were creat ed with the HealPix mapping 
tool (Gorski et ah 2005), which divides the sphere into 
pixels ot equal areas. For the present work, a pixelization 
parameter of Ngide = Id was used, which corresponds to 
3072 pixels in total, with a mean spacing of 3.67°. In 
Figure the excess regions, with respect to the average 
isotropic flux, are identified by a red color and deficit re¬ 
gions by a blue one. Therefore, a pixel in which many 
particles pass through will be represented in a stronger 
red color than one that has only a few events. 

Figure under the assumption of the Earth’s loca¬ 
tion in the interarm zone with Linj=100 pc, shows the 
sky maps progression of the 30 PeV cosmic ray arrival 
direction distribution with the dipole density gradient 
weight at different propagation distances of R = 10 pc, 
20 pc, 60 pc, and 90 pc. On each map, a dipole fit was 
performed, and the resulting dipole component was sub¬ 
tracted. Such a component may be different for each 
map, since it depends on the actual magnetic field struc¬ 
ture at a given distance. The sky maps show that by 
increasing the propagation distance the arrival direction 
distribution progressively develops smaller structures up 
the mean free path (60 pc in this case, see Table Q. At 
larger distances, the arrival direction distribution reaches 
a statistically steady configuration (in Figure]^ only the 
90 pc propagation distance is shown). Only propagation 
processes within the mean free path are responsible for 
the actual arrival directio n distribution of cosmi c rays at 
Earth, as discussed also in Giacinti & Sigl (2012). What¬ 
ever happens at larger distances is reshutiied by the scat¬ 
tering processes and is only relevant in the generation of 
the seed particle density gradient at large scale. 

The actual distribution of cosmic rays depends on the 
specific realization of the turbulent magnetic field and on 
the particle energy. Eigureffl shows the steady-state ar¬ 
rival direction distribution of750 TeV (on the left) and 30 
PeV (on the right) cosmic rays. This figure qualitatively 
shows that in the higher energy case the anisotropy of the 
distribution shows more small-scale angular regions than 
in the lower energy case. Eor a quantification of such a 
visual property, it is necessary to calculate the angular 
power spectrum. The map from 7.5 PeV presents almost 
the same distribution as the one for 750 TeV, since the 
only differences between them are the assumption on the 
injection scale and that they are independent sets (as 
described in Section 2.2 and shown in Table [^. The par¬ 


ticles at 7.5 PeV energy and Li^j = lOOpc are physically 
equivalent to particles at 750 TeV with Linj = lOpc. 

4.3. Angular Power Spectrum 

The sky maps of arrival direction distributions de¬ 
scribed in the previous sections are not dissimilar to 
experimental observations. However, it is the deter¬ 
mination of the angular power spectrum that makes it 
possible to quantify the formation of small-scale struc¬ 
ture anisotropy arising from scattering off magnetic tur¬ 
bulence. With the power spectrum, a sky map of ar¬ 
rival direction distribution is decomposed into spherical 
harmonics to provide information on the angular scale 
contribution of the anisotropy. The spectrum quan¬ 
titatively indicates which multipole moments £ in the 
spherical harmonic expansion contribute to the observed 
map. The IceGube observatory provided a power spec- 


TeV scale ( 

Abbasi et al. 

2011 Aartsen et al. 

the HAWC 

gamma-ray o 

bservatory provided 


Northern Hemisphere in the TeV scale (Abeysekara et al 


2014). The angular power spectrum is determined using 


the anafast tool in the HealPix framework. 

Eigure shows the angular power spectrum from the 
numerical trajectory integration set of 30 PeV protons 
at propagation distances from the initial unperturbed 
dipole anisotropy at 10 pc, 20 pc, 60 pc and 90 pc, as 
in Eigure ^ but without subtracting the dipole compo¬ 
nent. In the figure, the result from the observations made 
by the IceGube observatory is included as well. Note 
that the higher power values observed at low £, not re¬ 
produced in the calculations, are most probably due to 
the partial sky view of the IceGube observatory, which 
causes correlations with the largest scale spherical har¬ 
monic moments. The experimental observation is for a 
median cosmic ray energy of about 20 TeV, much lower 
than the numerical calculations. The numerical calcula¬ 
tion sets are normalized to the dipole component (i.e., 

= 1) of experimental power spectrum for the farthest 
propagation distance of 90 pc only. At shorter propa¬ 
gation distances, the normalization corresponds to the 
relative power obtained from the calculations. This nor¬ 
malization is valid since we are interested in the small- 
scale features, not on the assumptions on the large-scale 
anisotropy. 

Note that the numerical calculation shows, within sta- 
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Fig. 5. — Angular power spectrum of the arrival direction distribution of 30 PeV trajectories set of Tableland Figurewith dipole 
weight injected at a distance of 10 pc (in cyan), 20 pc (in green), 60 pc (in blue and corresponding to the mean free path) and of 90 pc 
(in red). The gray bands show the Icr, 2cr and 3cr b ands for a large set of isotropic sky maps. The black circles are the results from the 
IceCube observatory at a median energy of 20 TeV ([Santander et al.|2013|). Note the difference in energy scale between the experimental 
data and the numerical calculations. 




Fig. 6. — Angular power spectrum of the arrival direction distribution of 750 TeV (blue line on the left) and 30 PeV (blue line on the 
right) trajectori es sets of Tab le with dipole weight injected at the corresponding mean free path distance. The red line is the power 
spectrum from (|Ahlers||2014|). Tme black circles are the results from the Ice Cube observatory at a median energy of 20 TeV. The gray 
crosses and error bars sJiow the Icr band for a large set of isotropic sky maps ([Santander et al |2013|). Note the difference in energy scale 
between the experimental data and the numerical calculations. 


tistical fluctuations, that the power of the dipole compo¬ 
nent decreases with increasing propagation distance, due 
to the transfer of part of it into the higher I components. 
At very short distances, i.e., smaller than the mean free 
path, the low multipole moments are dominant. How¬ 
ever, as the distance increases, more power is transferred 
to the higher multipole moments, which is the signature 
of particle interactions with the turbulent magnetic field. 


Once they reach the mean free path distance, in this case 
60 pc, the power spectrum converges to a steady configu¬ 
ration, as the sky maps in Figure^ show. In Figurethe 
power spectrum of an isotropic nux of the same number 
of particles is included as a gray band and properly nor¬ 
malized Q It is evident that the distribution at distance 

^ The power spectrum of the isotropic ffux is calculated by gener- 
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Xmfp and beyond is not only a random distribution but 
also possesses a definite structure. It is noted that the an¬ 
gular power spectra calculated for propagation distance 
longer that 90 pc show the same relative normalization 
and shape of that corresponding to the mean free path. 
This is compatible with the achievement of steady state 
in the anisotropy structures as evident in the sky maps 
of Figure The approach of the power spectra to the 
band of isotropic distribution for > 30 means that, for 
the numerical realization studied here, the smaller angu¬ 
lar scales are indistinguishable from random fluctuations 
of the isotropic flux. 

As observed in Figure the angular structure gener¬ 
ated from the effect of scattering off magnetic turbulence 
depends on the particle energy (750 TeV and 30 PeV 
shown). This is evident also in the corresponding power 
angular spectra in Figure calculated at the propaga¬ 
tion distance of their mean free paths. The figure also 
shows the experimental result from the IceCube obser¬ 
vatory at 20 TeV median energy, with the corresponding 
power of the isotropic distribution background and with 
the curve expected fro m the hierarch ical breakup of an¬ 
gular components from Ahlers (|2014), normalized to the 
dipole component of the iceCiibe result. Note that the 
angular power spectrum at higher energy is flatter than 
that at lower energy. This is compatible with the exis¬ 
tence of more small-scale structures as evident in the sky 
maps of Figure 

5. DISCUSSION 

We have shown how small-scale anisotropy arises from 
the interaction of particles with the turbulent magnetic 
field. Specifically, we have addressed how the integration 
of trajectories in an MHD turbulent magnetic field pro¬ 
vides a realistic understanding of the small-scale features 
present in the observations of cosmic ray anisotropy at 
Earth. The ISM is in a plasma state, where the MHD 
equations serve as a model for its dynamics, therefore an 
MHD turbulent magnetic field is a natural approach to 
study the magnetic field prop erties in the ISM. Previous 
work (Giacinti & Sigl 2012) had considered the effects 
of magnetic turbulence on cosmic ray distribution. In 
this study, the authors have considered synthetic turbu¬ 
lence, which, on one hand, lacks the proper development 
of the gas dynamics but, on the other hand, provides the 
first qualitative connection between a magnetic turbulent 
field and cosmic ray arrival directions. In compressible 
MHD turbulence, scattering efficiency strongly depends 
on the wave type and how the particle gyroradius com¬ 
pares to the turbulence scales. Specifically, fast modes 
are identified as the main source for cosmic ray scatter¬ 
ing (|Yan & Lazarian||2008). 

The dynamics of the different turbulence modes and 
the relationship between particle energy and turbulence 
scale determine the actual scattering efficiency, which is 
most responsible for the breakout of a global cosmic ray 
density gradient into small-scale angular anisotropy. The 
angular power spectrum calculated for the data sets stud¬ 
ied in this work appears to evolve in time until particles 
cross the mean free path distance. In a steady state, the 


ating 10,000 realizations of uniform particle direction distribution, 
matching the number of the integrated particle trajectories, and 
calculating the power spectrum for each of them. The 68%, 90%, 
and 99% containment bands are reported in Figure 


shape of the angular power spectrum is a function of the 
magnetic field structure and of the consequent effects on 
particle propagation. 

Studying particle trajectories in MHD magnetic turbu¬ 
lence provides a more realistic framework to investigate 
the behavior of cosmic ray propagation in the ISM, where 
turbulence is expected to be anisotropic, although MHD 
turbulence simulations typically represent a significantly 
more restrictive inertial range than the actual astrophys- 
ical plasmas. 

The ISM is a complex environment and its exact rep¬ 
resentation is difficult to achieve; therefore, our MHD 
magnetic field can be considered one possible configura¬ 
tion of the magnetic field in the ISM. For this reason, 
direct comparison should be on the angular power spec¬ 
trum itself, not on the exact topology of the small-scale 
features in our maps. 

The framework for the application of Liouville’s the¬ 
orem is provided in the context of cosmic ray arrival 
distribution and applied through the back-propagation 
method. Although particle trajectories appear to suffer 
from mild deviation from adiabaticity, Liouville’s the¬ 
orem was applied in this particular case to study the 
first order effects of magnetic turbulence on the global 
topology of particle trajectories. This is because parti¬ 
cles do not experience severe scattering in their trajecto¬ 
ries, as shown in the first adiabatic invariant calculations; 
nonetheless, if the magnetic field were to vary dramat¬ 
ically with respect to the gyroradius of the particles, it 
prohibit application of Liouville’s theorem. The spread 
in the magnetic moment /i distribution in Figure [^sug¬ 
gests that some trajectories may have experienced more 
scattering than average in their propagation. This effect 
will manifest in the anisotropy through higher power at 
high since these particles will have had greater inter¬ 
actions with the turbulent field, resulting in a slightly 
flatter angular power spectrum. If the distribution on 
the (Jf^/jH plot had peaked at a higher value or a progres¬ 
sive trend had appeared on the mean magnetic moment 
plot of the data set, this would be an indication of a clear 
violation of Liouville’s theorem. Neither of these trends 
have statistically or significantly occurred in our calcu¬ 
lations; nevertheless, we would expect them to appear 
in magnetic fields that have a strong variation in scale 
on the order of the gyroradius of the particles. Future 
studies will need to use MHD turbulence simulation with 
wider inertial range to include enlarge the energy range 
in which magnetic turbulence affects cosmic ray distribu¬ 
tion. In this way, these studies will reveal the connection 
between the angular power spectrum of the cosmic ray 
arrival direction distribution and the turbulence proper¬ 
ties, naturally accounting for spurious effects derived by 
the numerical methods used. 

As mentioned in Section |2.2[ the trajectory back- 
propagation method is intended to provide a high effi¬ 
ciency in the studies of particle propagation in magnetic 
fields. Such a numerical method provides a mapping 
of the regions in space that are magnetically connected 
to the arrival point, where particles are assumed to be 
isotropically distributed. Therefore, appealing to the 
conservation of the total power across all spherical har¬ 
monic contributions, as dictated by Liouville’s theorem, 
makes higher multipole moments possible. Anisotropy 
was studied as a function of an initial large-scale density 
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gradient at a large distance from the arrival point. Such 
a density gradient, however, is generated by the same 
propagation processes that produce all angular struc¬ 
tures in the arrival direction distribution as well. In fact, 
smaller scale structures can arise in trajectory forward- 
propagation integration methods, where particle escape 
is naturally accounted for, even wi thout assuming an ini - 
tial global anisotropy, as shown in|Rettig & Pohl| (|2015). 
On the other hand, as long as a global anisotropy is devel- 
oped at some distance larger than Xmfp^ independently 
on how it is generated, the arrival direction distribution 
reaches a steady state angular power spectrum, and the 
effects of seed anisotropy and observed anisotropy can 
be disentangled. 

For the mean free path calculation, it is shown to be 
dependent on energy. In this regime, from 750 TeV to 30 
PeV, the Xmfp increases with energy. This is in agree¬ 
ment with how the more energetic particles interact with 
the perturbations of the magnetic field. In the case of 
the 750 TeV and 7.5 PeV sets, the minimum scale in the 
power spectrum of turbulence is of the size of the gyrora- 
dius of the particles, which may cause an overestimation 
of the Xmfp, since there is limited power in the physical 
perturbations that interact with the particles. For our 
30 PeV particles, we are well above this minimum scale, 
so the Xmfp is unaffected. 

The anisotropy maps and their corresponding angular 
power spectrum can tell us about the propagation of the 
particles themselves in the turbulent field. At very short 
distances, the low multipole moments are completely 
dominant and hold most of the power, as shown with 
the line of 10 pc in Figure The reason for this is that 
the particles have not had enough time to interact with 
the features of the magnetic field, and the distribution 
is still reminiscent of the initial configuration. However, 
the particles continue to interact and structures start to 
develop. As we can see with the line of 20 pc in Figure 
the higher multipole moments start to rise, with small- 
scale features becoming highly visible even in the maps 
in Figure One interesting feature is that the small- 
scale structures develop within the mean free path, but 
once they reache this scale, the maps do not change sig¬ 
nificantly. This observation is confirmed in the angular 
power spectrum. Therefore, the distribution of power be¬ 
tween the different multipole moments reaches a steady 
state. Of course, the particles keep moving and interact¬ 
ing after they have reached one Xmfp, but the anisotropy 
itself does not change. From Figure it is possible to 
see that this steady distribution is not isotropic, yet it 
possesses a definite structure that is dependent on the 
nature of the turbulent magnetic field. Consequently, 
the observed anisotropy is for the most part created in 
the last Xmfp of the particle’s trajectory, and it becomes 
a way to indirectly probe the local ISM. 

In Figure we can see from the comparison with the 
observations that the experimental data behaves accord¬ 
ing to what we would have expected from a lower en¬ 
ergy. In the case of 20 TeV, the distribution of power 
among the higher multipole moments is lower than at the 
higher energy, i.e., 30 PeV. One of the causes is that the 
30 PeV particles interact with perturbations that carry 
more power than the ones at a lower energy. There¬ 
fore, when an interaction process occurs with these per¬ 
turbations, the higher energy particles are affected more 


strongly and more small-scale structure is created. 

This effect is confirmed in Figure where the spec¬ 
trum at the Xmfp of energies 750 TeV and 30 PeV is 
compared with the observations at the median energy of 
20 TeV. The distribution at 750 TeV is similar to the ob¬ 
servations at 20 TeV; however, for 30 PeV, the distribu¬ 
tion of high multipole moments is much higher and flatter 
than that at lower energy. Therefore, greater small-scale 
anisotropy is observed at higher energies. This is also 
easily seen in the maps in Figure 

The flatter angular power spectrum obtained with 30 
PeV compared to 750 TeV protons tells us about dif¬ 
ferences in pitch angle scattering as a function of energy. 
However it is important to note that the 750 TeV data set 
corresponds to a gyroradius scale close to the dissipation 
region of the MHD turbulence inertial range. It is likely 
that scattering is underestimated, thus preventing a full 
development of small-scale angular structures in the par¬ 
ticle anisotropy. From this point of view, it is reasonable 
to assume that the resemblance of the 750 TeV proton 
power spectrum to experimental data should be consid¬ 
ered coincidental but still in agreement with the fact that 
lower energies should have less structure, as mentioned 
above. On the other hand, the difference of the 30 PeV 
proton power spectrum with experimental data does not 
necessarily mean that the fundamental processes respon¬ 
sible for the small-scale anisotropy are overestimated in 
the present study. The results show an energy depen¬ 
dence in the shape of the angular power spectrum that 
needs further study, requiring MHD turbulence simula¬ 
tion with a significantly wider inertial range. 

Figure includes the angular power spectrum result¬ 
ing from hierarchical decay of angul ar scales if Lion ville’s 


theorem is satisfied, as calculated by Ahlers (2014). Here 
the different curves have been norrnalized to the dip ole 
of the observational data, as discussed in Section 
that o nly the small-scale structure is relevant. In 


4.3 

, so 

Ah; 

ers| 

ves 

to 


(2014), an existing global dipole anisotropy evol 
create higher multipole moments. The difference be¬ 
tween that result and the one shown in the present work 
resides in the fact that here the shape of the distribution 
is determined by the specific turbulence characteristics 
of the magnetic field and the gyrorad ius of the part icles. 
The 750 TeV case appears similar to Ahlers (2014), but 
since this set is at the damping scale and scattering is 
likely to be underestimated, we were already expecting 
le ss structure to be present . 

Ahlers & Mertschl (2015) studied the effect of relative 


diftusion, i.e., from correlated nearby trajectories, as the 
contribution to the development of small scale anisotropy 
structures. In this complementary approach, the angular 
power spectrum is calculated based on the effect that a 
particle density gradient has on the trajectory topology 
shaped by homogeneous isotropic magnetic turbulence. 
In the present work, the relationship between the shape 
of angular power spectrum and the scattering processes 
induced by Alfvenic, slow, and fast modes in a compress¬ 
ible MHD turbulence were studied. This will make it 
possible to identify physical properties that shape the 
angular power spectrum, so that the problem can be in¬ 
verted, i.e., the angular power spectrum and different 
cosmic ray energies and masses can be used to probe the 
properties of the interstellar magnetic turbulence. 
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6. CONCLUSIONS 

This work explores the possibility that the local inter¬ 
stellar magnetic field could shape the arrival direction 
distribution of high-energy cosmic rays. An MHD tur¬ 
bulent magnetic field was used as a propagating medium 
that resembles the ISM, and particle trajectories with 
various energies were integrated in this field. To obtain 
the anisotropy in arrival direction distribution at Earth, 
the theoretical framework for application of Liouville’s 
theorem was provided, and the theorem was shown to be 
applicable in this specific case. Nonetheless, there could 
be cases where the magnetic field varies abruptly, which 
would ultimately violate the application of the theorem. 

The results presented in this work show that small- 
scale anisotropy arises from the interaction of cosmic rays 
with the local turbulent magnetic field, as confirmed in 
the sky maps and angular power spectra. The angular 
power spectrum becomes flatter the higher the energy; 
therefore, experimental data at the 10 TeV scale is ex¬ 
pected to be steeper than the numerical calculations pre¬ 
sented in this paper. Cosmic rays with less rigidity are 
more sensitive to the lower power small-scale magnetic 
perturbations. Since cosmic rays are not composed only 
of protons, the contribution of heavy ions would yield a 
steeper angular power spectrum. 

The inertial range of our turbulent magnetic field pro¬ 
vides a limitation on the lowest energy that could be 
studied, in this case 750 TeV. Therefore, in the future, it 
will be necessary to extend this inertial range and sam¬ 
ple even lower energies, so that direct comparisons with 


the observations at 20 TeV will be possible. Still, it is 
expected that at these TeV energies the effects of the he¬ 
liosphere should be present as well, which would provide 
a primary topic for future study and publication. An¬ 
other factor to improve is the number of particles that 
are propagated. If we were to have at least 10^ events, 
it would be possible to obtain a clearer view of the dis¬ 
tribution of cosmic rays at Earth. Investigating how the 
properties of our local turbulent magnetic field might in¬ 
fluence TeV-PeV cosmic ray arrival direction distribution 
will provide the basis for further exploring the observed 
anisotropy, and it will open the doors for a better under¬ 
standing of our local interstellar medium. 

Many thanks to Markus Ahlers, Juan Carlos Dfaz 
Velez, and other colleagues at WIPAC, and to Ellen 
Zweibel, Eederico Eraschetti, Martin Pohl, Rahul Kumar 
and Jungyeon Cho for useful and constructive discus¬ 
sions on the topic. AL acknowledges the support of NSE 
grant AST 0808118, NASA grant X5166204101, and the 
NSE-sponsored Center for Magnetic Self-Organization. 
PD acknowledges support from WIPAC and the U.S. 
National Science Eoundation-Office of Polar Programs. 
VLB thanks the support of the Brazil-U.S. Physics Ph.D. 
Student and Post-doc Visitation Program. This work 
was partially supported by the Research Experiences for 
Undergraduates (REU) Program of the National Science 
Eoundation under Award Number AST-1004881. 


APPENDIX 


NUMERICAL APPROACH AND ACCURACY 

Parti cle t rajectories are calculated by integrating the 6D set of equations of motion, eqs. 1 and As stated in 
section the integration is performed numerically using the Bulirsch-Stoer method, which is co nsidered one of the 
best known integratio n algorithms satisfying both high accuracy and effi ciency (Press et al.||1986) and widely applied 
in the literature (e.g., jGiacalone & Jokipii| ( 1999|) andjXu & Yan| (|2013|)). The Bulirsch-Stoer integration algorithm, 
is a known method for numerical calculation of ordinary difterenUal equation solutions, that combines the so-called 
Richardson extrapolation (to improve the rate of convergence of a sequence) and the modified midpoint method 
(which advances a vector of dependent variables y(x) from a point x to a point x -L H by a sequence of n substeps 
each of size h). The result is that Bulirsch-Stoer algorithm provides high accuracy with relatively low computational 
effort. The accuracy of the algorithm is further controlled, during the numerical calculation, by monitoring the local 
truncation error estimated at each time step. If the relative error is larger than the relative tolerance level of 10“^, the 
step size is adaptively reduced in order to limit the error accumulation in both momentum and spatial coordinates, 
across the maximum integration time used in this work (corresponding to no more than 10,000 gyrations). Such error 
accumulation needs to be monitored for each specific problem in which the integration algorithm is used. In particular, 
the accuracy on the spatial and momentum coordinates were studied. 

The performance of the Bulirsch-Stoer algorithm on the accuracy of momentum coordinates in our numerical cal¬ 
culation can be seen in Eigurej^ where the particle energy variation (due to loss of accuracy) from that at t = 0 is 
plotted as a function of the nurnber of gyro-orbits Uq t (for a single particle and for the average of all particles used in 
our 30 PeV sample). 

The maximum relative mean deviation from perfect energy conservation is found to be about 8.5x10“^ for the 
sample used in this work (100,000 particles in the 30 PeV energy set), although a single particle can reach a violation 
at the 1.6x10“^ level. This precision level in the energy conservation guarantees that particle trajectories are not 
significantly affected by numerical accuracy limitations, which are found to be marginal under the conditions of this 
study. 

Since the adaptive time step algorithm constrains both spatial and momentum coordinates to the same relative 
error level, the accuracy in spatial coordinates is vl, even after 10,000 gyrations. Numerical diffusion, therefore, is 
limited to a level much smaller than Bohm diffusion, which i s several orders of mag nitude below diffusion induced by 
the stochastic wandering of magnetic field lines at all scales (Laza r ian fc Yan 2014). 

(|2013[ ), which used the same integration stepping 


The effect on the particle set size was assessed in Xu V Yan 


algorithms as in this work. In that paper, the perpendicular diffusion coefficient becomes stable when the sample 
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Fig. 7. — Accuracy of the conservation of energy for a single particle (on the left) and for the average particle sample (on the right) of 
the 30 PeV set of Table [T] 

size reaches about 1000 particles. The sets used in this study contain 100,000 or more particles (see Table Q, thus 
minimizing statistical accuracy effects on the global behavior of the ensemble of particles. 
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